Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany
Bashar Jarayseh
Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany
Ailsa Howard
South Bay Banded Dotterel Project, Kaikoura, New Zealand
Ted Howard
South Bay Banded Dotterel Project, Kaikoura, New Zealand
Tony Habraken
Port Waikato Banded Dotterel Project, Port Waikato, New Zealand
Emma Williams
Department of Conservation, Christchurch, New Zealand
Colin O`Donnell
Department of Conservation, Christchurch, New Zealand
Bart Kempenaers
Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany
Published
January 27, 2025
Main conclusions since last meeting in Nov, 2024
Sex difference in maximum rufous band score (strong effect)
Males have a more full and immaculate rufous breast band than females
No age effect on maximum rufous band score
Migratory-status difference in maximum rufous band score (weak, albiet significant effect)
Residents of both sexes have slightly more full and immaculate rufous breast bands than migrants
No sex-specific timing of the pre-basic moult
No migratory-status-specific timing of the pre-basic moult
Age-specific timing of pre-basic moult (strong effect)
younger birds initiate their pre-basic moult earlier than older birds
No effect of maximum rufous band score on prolonged breeding
No effect of maximum rufous band score on early breeding
No effect of prolonged breeding on maximum rufous band score
Positive relationship between the date when an individual was last seen with its maximum breeding plumage and the date when it was last observed breeding
No age effect on prolonged breeding
No age effect on early breeding
Exploration of moult data extracted from photos
Note: data last modified 7-Nov-2024
This dataset contains the scored moult data for all usable photos from the 2021-2022, 2022-2023, and 2023-2024 breeding seasons at Kaikoura
summarize the number of individuals in the dataset
[1] 94
assess how many individuals have more than a single season of data
[1] 63
summarise the number of usable photos for each individual across the three seasons
the final dataset
Code
dat %>%datatable(class ='cell-border stripe', rownames =FALSE, filter ='top')
Plot the sampling distribution of Ailsa’s photos across the seasons
The “score” value is based on a scale from 1 to 7 describing the immaculate-ness and full-ness of the rufous breast band seen in the photos (Kok, Hogan, and Piersma (2020))
NOTE FOR BASHAR (26-01-2025):
would be great to add a figure here showing example photos of an individual male and female over several phases of the moult. For example, two-columns with male on left, female on right, and 5 rows of images starting with their fullest and most immaculate state and ending with their basic plumage. Would be useful to show the scores you gave for these individuals in the figure. Based on the data, it seems like Male "RRGO" and Female "RWBB" would be a good candidates given that they recieved a wide range of scores
# A tibble: 197 × 4
# Groups: sex, rings_comb [95]
sex rings_comb season n_unique_scores
<chr> <chr> <dbl> <int>
1 M RRGO 2022 7
2 F RWBB 2021 6
3 M RGGY 2023 6
4 M RRGO 2021 6
5 M RRGO 2023 6
6 M RRWW 2021 6
7 F RGGL 2021 5
8 F RROB 2022 5
9 F RWRR 2023 5
10 M RROB 2022 5
# ℹ 187 more rows
2021-2022 season data
2022-2023 season data
2023-2024 season data
Sources of individual variation in the maximum score of the rufous band
Sex-specific differences in breeding plumage
data wrangle for sex model - note that the max plumage score recorded during the core breeding months (Aug-Dec) is extracted for each individual for each year (“max_breeding_score”).
Code
# extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual.ind_breeding_scores <- dat %>%mutate(score =as.numeric(score)) %>%mutate(breeding_season =ifelse(month(date) %in%c(8, 9, 10, 11, 12), 1, 0)) %>%filter(breeding_season ==1) %>%group_by(season, rings_comb, sex) %>%#filter(rings_comb == "RRBY") %>% select(score, date)summarise(max_breeding_score =max(score),date_of_max_score = date[which.max(score)],mode_breeding_score =get_mode(score)) %>%ungroup() %>%mutate(max_breeding_score_mean =round((max_breeding_score + mode_breeding_score)/2)) %>%mutate(date_J_of_max_score =as.numeric(format(date_of_max_score +181, "%j")))
Assess sample sizes of each sex
# A tibble: 2 × 2
sex `n_distinct(rings_comb)`
<chr> <int>
1 F 49
2 M 35
mixed model with individual and season as random intercepts and sex as fixed effect
Code
# linear mixed model for the difference in max score between sexesmod_max_score <-lmer(max_breeding_score ~ sex + (1| rings_comb) + (1| season), data = ind_breeding_scores)
Males tend to have a more full and immaculate rufous bands than females, however there is substantial overlap in their distributions.
sex-specific breeding plumage distributions
effect-size table for sex-specific breeding plumage model
Banded Dotterel rufous band score
Mean estimate
95% confidence interval
Fixed effects 𝛽 (standardized)
Male (Sex)
1.01
[0.73, 1.27]
Partitioned 𝑹²
Total Marginal 𝑹²
0.29
[0.17, 0.41]
Sex
0.29
[0.17, 0.41]
Total Conditional 𝑹²
0.51
[0.33, 0.66]
Random effects 𝜎²
Individual
0.17
[0.04, 0.31]
Year
0.03
[0, 0.13]
Residual
0.42
[0.31, 0.57]
Adjusted repeatability 𝑟
Individual
0.27
[0.07, 0.46]
Year
0.04
[0, 0.17]
Residual
0.69
[0.49, 0.9]
Sample sizes 𝑛
Years
3
Individuals
83
Observations
164
sex-specific breeding plumage model predictions
Age- and sex-specific variation in breeding plumage (of known-aged birds)
data wrangle for sex-age model. Standardize max_breeding_score according to sex to make males and females to make them comparable across ages (i.e., does max_breeding_score change across age while accounting for underlying sex differences?)
Code
# extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual.ind_breeding_scores_age <- dat %>%mutate(score =as.numeric(score)) %>%mutate(breeding_season =ifelse(month(date) %in%c(8, 9, 10, 11, 12), 1, 0)) %>%filter(breeding_season ==1) %>%group_by(season, rings_comb, sex, age) %>%summarise(max_breeding_score =max(score)) %>%ungroup() %>%group_by(sex) %>%mutate(std_max_breeding_score = (max_breeding_score -mean(max_breeding_score, na.rm =TRUE)) /sd(max_breeding_score, na.rm =TRUE)) %>%filter(!is.na(age))
Assess sample sizes of each sex for individuals with known age
# A tibble: 2 × 2
sex `n_distinct(rings_comb)`
<chr> <int>
1 F 14
2 M 10
mixed model of standardized max breeding score (by sex) with individual as random intercept and age as fixed effect (Note: we tried to include season as an additional random intercept, but the low sample size of repeated measures across seasons caused convergence issues of the model)
Code
# linear mixed model for the difference in max score between sexes across agemod_max_score_age <-lmer(std_max_breeding_score ~ age + (1| rings_comb), data = ind_breeding_scores_age)
effect-size table for sex- and age-specific variation in breeding plumage
no effect of age on the maximum extent of an individual’s rufous band score while controlling for underlying sex-specific variation
Banded Dotterel rufous band score
Mean estimate
95% confidence interval
Fixed effects 𝛽 (standardized)
Age
0.09
[-0.24, 0.37]
Partitioned 𝑹²
Total Marginal 𝑹²
0.01
[0, 0.13]
Age
0.01
[0, 0.13]
Total Conditional 𝑹²
0.38
[0.01, 0.72]
Random effects 𝜎²
Individual
0.27
[0, 0.66]
Residual
0.45
[0.2, 0.75]
Adjusted repeatability 𝑟
Individual
0.37
[0, 0.72]
Residual
0.63
[0.28, 1]
Sample sizes 𝑛
Years
3
Individuals
24
Observations
42
age-specific breeding plumage model predictions
Migratory status- and sex-specific variation in breeding plumage (of known-status birds)
data wrangle for sex-migratory status model
Code
# extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual.ind_breeding_scores_status <- dat %>%mutate(score =as.numeric(score)) %>%mutate(breeding_season =ifelse(month(date) %in%c(8, 9, 10, 11, 12), 1, 0)) %>%filter(breeding_season ==1) %>%group_by(season, rings_comb, sex, age, migratory_status) %>%summarise(max_breeding_score =max(score)) %>%ungroup() %>%group_by(sex) %>%mutate(std_max_breeding_score = (max_breeding_score -mean(max_breeding_score, na.rm =TRUE)) /sd(max_breeding_score, na.rm =TRUE)) %>%filter(migratory_status %in%c("M", "R"))
Assess sample sizes of each sex and migratory status (for birds with known status)
Code
# Assess sample sizes of each sex and season # (13 females (2 migrant, 11 resident), 8 males (3 migrant, 5 resident))ind_breeding_scores_status %>%group_by(sex, migratory_status) %>%summarise(n_distinct(rings_comb))
# A tibble: 4 × 3
# Groups: sex [2]
sex migratory_status `n_distinct(rings_comb)`
<chr> <chr> <int>
1 F M 2
2 F R 11
3 M M 3
4 M R 5
mixed model of standardized max breeding score (by sex) with individual as random intercept and migratory status as fixed effect (Note: we tried to include season as an additional random intercept, but the low sample size of repeated measures across seasons caused convergence issues of the model)
Code
# linear mixed model for the difference in max score between sexesmod_max_score_status <-lmer(std_max_breeding_score ~ migratory_status + (1| rings_comb), data = ind_breeding_scores_status)
effect-size table and plot for sex- and migratory-status breeding plumage
no effect of migratory status on an individual’s rufous band score while controlling for underlying sex-specific variation
Banded Dotterel rufous band score
Mean estimate
95% confidence interval
Fixed effects 𝛽 (standardized)
Resident (Migratory Status)
0.55
[0.04, 1.16]
Partitioned 𝑹²
Total Marginal 𝑹²
0.10
[0, 0.34]
Migratory Status
0.10
[0, 0.34]
Total Conditional 𝑹²
0.23
[0.01, 0.61]
Random effects 𝜎²
Individual
0.08
[0, 0.38]
Residual
0.51
[0.23, 0.76]
Adjusted repeatability 𝑟
Individual
0.16
[0, 0.54]
Residual
0.84
[0.46, 1]
Sample sizes 𝑛
Years
3
Individuals
21
Observations
42
migratory status-specific breeding plumage model predictions
Sources of individual variation in the timing of pre-basic moult
Within-individual moult dynamics across sexes
wrangle data by calculating the individual proportional moult scores by comparing each score to a given individual’s max (determined in previous chunk). Piece apart sex-specific variation in the timing of pre-basic moult
Assess sample sizes of each sex
# A tibble: 2 × 2
sex `n_distinct(rings_comb)`
<chr> <int>
1 F 49
2 M 32
mixed effects binomial model assessing sex-differences in the timing of pre-basic moult while accounting for bird- and year-specific random variation in the timing of pre-basic moult (i.e., random slope model of date for individual and year).
Code
mod_moult_sex_bird_RIS <- lme4::glmer(prop_molt_score_max ~ date_J + sex + (date_J | ring_season),data = ind_prop_molt_scores, family = binomial)
Characteristic
log(OR)1
95% CI1
p-value
(Intercept)
22
18, 26
<0.001
Date
-0.10
-0.12, -0.08
<0.001
Sex
F
—
—
M
-0.03
-0.80, 0.75
>0.9
1 OR = Odds Ratio, CI = Confidence Interval
Within-individual moult dynamics across migratory status
# A tibble: 2 × 2
migratory_status `n_distinct(rings_comb)`
<chr> <int>
1 M 5
2 R 16
Code
# mixed effects binomial model comparing migratory status and date effect on the changes in moult scoresmod_moult_status_bird_RIS <- lme4::glmer(prop_molt_score_max ~ date_J + migratory_status + (date_J|ring_season),data = ind_prop_molt_scores_status, family = binomial)
# mixed effects binomial model comparing migratory status and date effect on the changes in moult scoresmod_moult_age_bird_RIS <- lme4::glmer(prop_molt_score_max ~ date_J + age + (date_J | ring_season),data = ind_prop_molt_scores_age, family = binomial)
Characteristic
log(OR)1
95% CI1
p-value
(Intercept)
53
18, 88
0.003
Date
-0.28
-0.47, -0.10
0.003
Age
2.2
0.61, 3.7
0.006
1 OR = Odds Ratio, CI = Confidence Interval
inspect the random slopes against the raw data from photos
Breeding status analysis
Relationships between breeding schedule and plumage
effect of breeding plumage score on prolonged breeding
do individuals with a fuller breeding plumage (for their sex) breed later than the rest of the population? No effect
ggplot() +geom_line(data =as.data.frame(allEffects(mod_last_breeding_date_prop_molt_mid))$last_breeding_date_J_snum,aes(x = last_breeding_date_J_snum, y = fit)) +geom_ribbon(data =as.data.frame(allEffects(mod_last_breeding_date_prop_molt_mid))$last_breeding_date_J_snum,aes(x = last_breeding_date_J_snum, ymax = upper, ymin = lower),lwd =1, alpha =0.25) +geom_point(data = ind_prop_molt_scores_last_breeding_date,aes(x = last_breeding_date_J_snum, y = mid_point_date_J),alpha =1, shape =20) +theme(text =element_text(family ="Verdana"),axis.ticks =element_blank(),plot.margin =unit(c(0.2,0.2,0.2,0.2), "cm"),panel.spacing =unit(1, "cm"),panel.background =element_rect(fill ='transparent', color =NA),plot.background =element_rect(fill ='transparent', color =NA), panel.grid.major.y =element_line(size =0.25, lineend ="round", colour ="grey90"),panel.grid.minor =element_blank(),panel.grid.major.x =element_line(size =0.25, lineend ="round", colour ="grey90")) +xlab("year-scaled date of last breeding attempt") +ylab("mid-point date between\nlast max score photo and first photo during moult")
effect of prolonged breeding on prolonged retention of full plumage
do individuals with late breeding attempts retain their full breeding plumage longer? There is a positive relationship between the date when an individual was last seen with its maximum breeding plumage and the date when it was last observed breeding.
ggplot() +geom_line(data =as.data.frame(allEffects(mod_age_first_breeding_date))$age,aes(x = age, y = fit)) +geom_ribbon(data =as.data.frame(allEffects(mod_age_first_breeding_date))$age,aes(x = age, ymax = upper, ymin = lower),lwd =1, alpha =0.25) +geom_point(data = ind_prop_molt_scores_first_breeding_date,aes(x = age, y = first_breeding_date_J_snum),alpha =1, shape =20) +theme(text =element_text(family ="Verdana"),axis.ticks =element_blank(),plot.margin =unit(c(0.2,0.2,0.2,0.2), "cm"),panel.spacing =unit(1, "cm"),panel.background =element_rect(fill ='transparent', color =NA),plot.background =element_rect(fill ='transparent', color =NA), panel.grid.major.y =element_line(size =0.25, lineend ="round", colour ="grey90"),panel.grid.minor =element_blank(),panel.grid.major.x =element_line(size =0.25, lineend ="round", colour ="grey90")) +xlab("age") +ylab("year-scaled date of first breeding attempt")
Meeting notes:
13-Nov-24: Bart, Luke, and Bashar via zoom
add the breeding data to the analysis: for each individual include the date of last breeding (i.e., last date seen with nest or brood) as a predictor of molt timing. Also do an analysis relating the timing of last breeding to the age of the individual. Bashar will send Luke the breeding data he received from Ted.
drop pre-October data from the molt timing analysis.
next steps: Bashar to eventually write up this study as a report (which could eventually be used as his MSc thesis if he wants). Priority will be on him writing up a proposal for his thesis (which has a concrete deadline and is graded). Bashar will present the study at our weekly Monday meeting in the first few months of 2025.
References
Kok, Eva MA, Jerry A Hogan, and Theunis Piersma. 2020. “Experimental Tests of a Seasonally Changing Visual Preference for Habitat in a Long-Distance Migratory Shorebird.”Ethology 126 (7): 681–93.
Source Code
---title: "Banded Dotterel Moult Study"subtitle: | Exploration of datasetdate: "`r format(Sys.time(), '%H:%M %d %B, %Y')`"author: - name: Luke Eberhart-Hertel orcid: 0000-0001-7311-6088 email: luke.eberhart@bi.mpg.de url: https://www.bi.mpg.de/person/115852/2867 affiliations: - ref: bk - name: Bashar Jarayseh affiliations: - ref: bk - name: Ailsa Howard affiliations: - ref: ah - name: Ted Howard affiliations: - ref: ah - name: Tony Habraken affiliations: - ref: th - name: Emma Williams affiliations: - ref: ew - name: Colin O`Donnell affiliations: - ref: ew - name: Bart Kempenaers affiliations: - ref: bkaffiliations: - id: bk number: 1 name: Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany - id: ah number: 2 name: South Bay Banded Dotterel Project, Kaikoura, New Zealand - id: th number: 2 name: Port Waikato Banded Dotterel Project, Port Waikato, New Zealand - id: ew number: 3 name: Department of Conservation, Christchurch, New Zealandformat: html: toc: true code-fold: true code-tools: true self-contained: true highlight-style: github theme: Cosmoexecute: warning: false cache: trueeditor_options: chunk_output_type: consolebibliography: references.bib---```{r, include=FALSE}knitr::opts_chunk$set(cache = TRUE)``````{r, message=FALSE, results='hide', warning=FALSE, cache=FALSE, include=FALSE}## Prerequisites### R packages# - The following packages are needed for analysis and can be easily installed from [CRAN](http://cran.r-project.org/) or GitHub by running the following code chunk:# a vector of all the packages needed in the projectpackages_required_in_project <- c("tidyverse", "readxl", "RMark", "RColorBrewer", "patchwork", "mapview", "lubridate", "extrafont", "here", "DT", "leaflet", "sf", "leafpop", "tsibble", "corrplot", "gghalves", "gam", "pscl", "gamlss", "gt", "lme4", "ggpattern", "gtsummary", "effects", "lattice", "rptR", "partR2", "broom.mixed", "forcats", "glmmTMB", "sysfonts", "showtext", "patchwork", "standardize")# of the required packages, check if some need to be installednew.packages <- packages_required_in_project[!(packages_required_in_project %in% installed.packages()[,"Package"])]# install all packages that are not locally availableif(length(new.packages)) install.packages(new.packages)# load all the packages into the current R sessionlapply(packages_required_in_project, require, character.only = TRUE)# set the home directory to where the project is locally based (i.e., to find # the relevant datasets to import, etc.here::set_here()``````{r, message=FALSE, results='hide', warning=FALSE, include=FALSE}### Plotting themes# - The following plotting themes, colors, and typefaces are used throughout the project:# Find fonts from computer that you want. Use regular expressions to do this# For example, load all fonts that are 'verdana' or 'Verdana'extrafont::font_import(pattern = "[V/v]erdana", prompt = FALSE) # check which fonts were loadedextrafont::fonts()extrafont::fonttable()extrafont::loadfonts() # load these into R# # Add Google font Verdana# font_add_google("Lato", "Lato")# # # Automatically use showtext to render text# showtext_auto()# define the plotting theme to be used in subsequent ggplotsluke_theme <- theme_bw() + theme( text = element_text(family = "Verdana"), legend.title = element_text(size = 10), legend.text = element_text(size = 8), axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 8), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 8), strip.text = element_text(size = 10), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_line(linewidth = 0.5, colour = "grey40"), axis.ticks.length = unit(0.2, "cm"), panel.border = element_rect(linetype = "solid", colour = "grey"), legend.position.inside = c(0.1, 0.9) )``````{r, include=FALSE}# Define a custom mode functionget_mode <- function(x) { uniq_vals <- unique(x) uniq_vals[which.max(tabulate(match(x, uniq_vals)))]}```## Main conclusions since last meeting in Nov, 2024#### Sex difference in maximum rufous band score (strong effect)Males have a more full and immaculate rufous breast band than females#### No age effect on maximum rufous band score#### Migratory-status difference in maximum rufous band score (weak, albiet significant effect)Residents of both sexes have slightly more full and immaculate rufous breast bands than migrants#### No sex-specific timing of the pre-basic moult#### No migratory-status-specific timing of the pre-basic moult#### Age-specific timing of pre-basic moult (strong effect)younger birds initiate their pre-basic moult earlier than older birds#### No effect of maximum rufous band score on prolonged breeding#### No effect of maximum rufous band score on early breeding#### No effect of prolonged breeding on maximum rufous band score#### Positive relationship between the date when an individual was last seen with its maximum breeding plumage and the date when it was last observed breeding#### No age effect on prolonged breeding#### No age effect on early breeding## Exploration of moult data extracted from photosNote: data last modified 7-Nov-2024This dataset contains the scored moult data for all usable photos from the 2021-2022, 2022-2023, and 2023-2024 breeding seasons at Kaikoura```{r, include = FALSE}# import data with all columns as character (so that no auto-formatting is done by R)dat <- read.csv("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Dataset_Molt_Final - modified 7.11.2024.csv", colClasses = "character") %>% mutate(Rings_comb = str_sub(Rings_comb, 1, 4)) %>% rename_with(tolower) # convert all column headers to lowercase# check values of Date column for mistakes...looks goodunique(dat$date)# check values of score column for mistakes...one row to fixunique(dat$score)# one row contains no data for the score (row 767)...needs to be fixeddat %>% filter(score == "")# check values of rings_comb column for mistakes...need to remove the non-unique combos (rows with XB and XW)unique(dat$rings_comb)dat %>% filter(rings_comb %in% c("XB", "XW"))dat <- dat %>% filter(!grepl("XB", rings_comb) & !grepl("XW", rings_comb))# import sex, migratory status, and age at banding informationdat_sexes <- read.csv("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/sex_status_banding.csv", colClasses = "character") %>% mutate(Rings_comb = str_sub(Rings_comb, 1, 4)) %>% rename_with(tolower) %>% # convert all column headers to lowercase distinct()# checkout the sex-type data# number of birds with sex-type datadat_sexes %>% pull(rings_comb) %>% unique() %>% length()# identify duplicates...none, gooddat_sexes[which(duplicated(dat_sexes)), ]# assess all birds with migratory status datadat_sexes %>% filter(migratory_status %in% c("R", "M"))# mutate the Date column into a date variabledat <- dat %>% mutate(date = paste(substring(dat$date, first = 7, last = 10), substring(date, first = 4, last = 5), substring(date, first = 1, last = 2), sep = "-") %>% as.Date()) %>% # subset to data with Molt == 1 filter(molt == 1) %>% # specify the season as the first calender year mutate(season = ifelse(month(date) < 7, year(date) - 1, year(date))) %>% # change to Julian date shifted for the Southern Hemisphere (1 = July 1) mutate(date_J = as.numeric(format(date + 181, "%j"))) %>% # join the sexes provided by Ailsa (note many-to-many is fine) left_join(., dat_sexes, by = "rings_comb", relationship = "many-to-many") %>% # remove individuals with unknown sex filter(sex != "U") %>% filter(sex != "N") %>% # add a ranking variable to sort the facets of the sampling distribution plots group_by(rings_comb) %>% mutate(n_photos = n(), n_scores = n_distinct(score)) %>% arrange(desc(n_scores)) %>% mutate( banding_date = paste(substring(banding_date, first = 7, last = 10), substring(banding_date, first = 4, last = 5), substring(banding_date, first = 1, last = 2), sep = "-") %>% as.Date(), season_ringed = ifelse(age_at_banding == "P" , ifelse(month(banding_date) < 7, year(banding_date) - 1, year(banding_date)), NA), age = ifelse(age_at_banding == "P", as.numeric(season - season_ringed), NA))# i should find a better way to calculate age```check max and min dates for each season```{r, echo= FALSE}# check max and min dates for each seasondat %>% group_by(season) %>% summarise(n_obs = n(), min_date = min(date), max_date = max(date))```summarize the number of individuals in the dataset```{r, echo = FALSE}dat %>% pull(rings_comb) %>% unique() %>% length()```assess how many individuals have more than a single season of data```{r, echo = FALSE}dat %>% select(rings_comb, season) %>% distinct() %>% group_by(rings_comb) %>% summarise(n_seasons = n()) %>% arrange(desc(n_seasons)) %>% filter(n_seasons > 1) %>% nrow()```summarise the number of usable photos for each individual across the three seasons```{r, echo=FALSE}dat %>% group_by(rings_comb, season) %>% summarise(n_photos = n()) %>% arrange(desc(n_photos)) %>% pivot_wider(., names_from = season, values_from = n_photos) %>% datatable(class = 'cell-border stripe', rownames = FALSE, filter = 'top')```the final dataset```{r, echo=TRUE}dat %>% datatable(class = 'cell-border stripe', rownames = FALSE, filter = 'top')```### Plot the sampling distribution of Ailsa's photos across the seasonsThe "score" value is based on a scale from 1 to 7 describing the immaculate-ness and full-ness of the rufous breast band seen in the photos (@kok2020experimental)### `NOTE FOR BASHAR (26-01-2025)`: `would be great to add a figure here showing example photos of an individual male and female over several phases of the moult. For example, two-columns with male on left, female on right, and 5 rows of images starting with their fullest and most immaculate state and ending with their basic plumage. Would be useful to show the scores you gave for these individuals in the figure. Based on the data, it seems like Male "RRGO" and Female "RWBB" would be a good candidates given that they recieved a wide range of scores````{r}dat %>%select(sex, rings_comb, season, score, date) %>%group_by(sex, rings_comb, season) %>%summarise(n_unique_scores =n_distinct(score)) %>%arrange(desc(n_unique_scores))```#### 2021-2022 season data```{r, fig.height=9, echo=FALSE}# 2021-2022 season dataggplot(data = dat %>% filter(season == 2021 & sex == "F") %>% mutate(score = as.numeric(score), rings_comb_ = reorder(rings_comb, n_photos, decreasing = TRUE))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ rings_comb_, ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) + xlab("week") + ggtitle("Females from the 2021-2022 breeding season in Kaikoura")``````{r, fig.height=8, echo=FALSE}ggplot(data = dat %>% filter(season == 2021 & sex == "M") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) + xlab("week") + ggtitle("Males from the 2021-2022 breeding season in Kaikoura")```#### 2022-2023 season data```{r, fig.height=10, echo=FALSE}# 2022-2023 season data (not as good coverage as the previous season)ggplot(data = dat %>% filter(season == 2022 & sex == "F") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) + xlab("week") + ggtitle("Females from the 2022-2023 breeding season in Kaikoura")``````{r, fig.height=7, echo=FALSE}ggplot(data = dat %>% filter(season == 2022 & sex == "M") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) + xlab("week") + ggtitle("Males from the 2022-2023 breeding season in Kaikoura")```#### 2023-2024 season data```{r, fig.height=10, echo=FALSE}# 2023-2024 season data ggplot(data = dat %>% filter(season == 2023 & sex == "F") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2023-07-01"), as.Date("2024-05-01"))) + xlab("week") + ggtitle("Females from the 2023-2024 breeding season in Kaikoura")``````{r, fig.height=6, echo=FALSE}ggplot(data = dat %>% filter(season == 2023 & sex == "M") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2023-07-05"), as.Date("2024-05-01"))) + xlab("week") + ggtitle("Males from the 2023-2024 breeding season in Kaikoura")```## Sources of individual variation in the maximum score of the rufous band### Sex-specific differences in breeding plumagedata wrangle for sex model - note that the max plumage score recorded during the core breeding months (Aug-Dec) is extracted for each individual for each year ("max_breeding_score").```{r}# extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual.ind_breeding_scores <- dat %>%mutate(score =as.numeric(score)) %>%mutate(breeding_season =ifelse(month(date) %in%c(8, 9, 10, 11, 12), 1, 0)) %>%filter(breeding_season ==1) %>%group_by(season, rings_comb, sex) %>%#filter(rings_comb == "RRBY") %>% select(score, date)summarise(max_breeding_score =max(score),date_of_max_score = date[which.max(score)],mode_breeding_score =get_mode(score)) %>%ungroup() %>%mutate(max_breeding_score_mean =round((max_breeding_score + mode_breeding_score)/2)) %>%mutate(date_J_of_max_score =as.numeric(format(date_of_max_score +181, "%j")))``````{r, include=FALSE}ind_breeding_scores %>% mutate(score_diff = max_breeding_score - mode_breeding_score) %>% ggplot() + geom_histogram(aes(score_diff))ind_breeding_scores %>% mutate(score_diff = max_breeding_score - mode_breeding_score) %>% arrange(desc(score_diff)) %>% select(season, rings_comb, sex, max_breeding_score, mode_breeding_score, date_of_max_score, score_diff)# three individuals with a difference in max and mode score greater than 2dat %>% filter(rings_comb == "RWBB" & season == "2021") %>% select(rings_comb, season, date, sex, n_photos, score) %>% arrange(date)dat %>% filter(rings_comb == "RWBG" & season == "2021") %>% select(rings_comb, season, date, sex, n_photos, score) %>% arrange(date)dat %>% filter(rings_comb == "RRWY" & season == "2022") %>% select(rings_comb, season, date, sex, n_photos, score) %>% arrange(date)dat %>% filter(rings_comb == "RRBY" & season == "2021") %>% select(rings_comb, season, date, sex, n_photos, score) %>% arrange(date)```Assess sample sizes of each sex```{r, echo=FALSE}ind_breeding_scores %>% group_by(sex) %>% summarise(n_distinct(rings_comb))```mixed model with individual and season as random intercepts and sex as fixed effect```{r, fold-show}# linear mixed model for the difference in max score between sexesmod_max_score <- lmer(max_breeding_score ~ sex + (1 | rings_comb) + (1 | season), data = ind_breeding_scores)``````{r, eval=FALSE, include=FALSE}tbl_regression(mod_max_score, intercept = TRUE, label = list(sex ~ "Sex"))# Derive confidence intervals of effect sizes from parametric bootstrappingtidy_mod_max_score <- tidy(mod_max_score, conf.int = TRUE, conf.method = "boot", nsim = 1000)# run rptR to obtain repeatabilities of random effectsrpt_mod_max_score <- rpt(max_breeding_score ~ sex + (1 | rings_comb) + (1 | season), grname = c("rings_comb", "season", "Fixed"), data = ind_breeding_scores, datatype = "Gaussian", nboot = 1000, npermut = 1000, ratio = TRUE, adjusted = TRUE, ncores = 4, parallel = TRUE)# run partR2 on each model to obtain marginal R2, parameter estimates, and beta# weightsR2m_mod_max_score <- partR2(mod_max_score, partvars = c("sex"), R2_type = "marginal", nboot = 1000, CI = 0.95, max_level = 1)R2c_mod_max_score <- partR2(mod_max_score, partvars = c("sex"), R2_type = "conditional", nboot = 1000, CI = 0.95, max_level = 1)stats_mod_max_score <- list(mod = mod_max_score, tidy = tidy_mod_max_score, rptR = rpt_mod_max_score, partR2m = R2m_mod_max_score, partR2c = R2c_mod_max_score, data = ind_breeding_scores)``````{r, include=FALSE}# saveRDS(stats_mod_max_score, file = "out/stats_mod_max_score.rds")stats_mod_max_score <- readRDS(file = "/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/out/stats_mod_max_score.rds")``````{r, include=FALSE}#### Table of effect sizes ----# Retrieve sample sizessample_sizes <- ind_breeding_scores %>% ungroup() %>% summarise(Year = n_distinct(season), Individual = n_distinct(rings_comb), Observations = nrow(.))sample_sizes <- as.data.frame(t(as.data.frame(sample_sizes))) %>% rownames_to_column("term") %>% rename(estimate = V1) %>% mutate(stat = "n")# clean model component namesmod_comp_names <- data.frame(comp_name = c("Male (Sex)", "Total Marginal \U1D479\U00B2", "Sex", "Total Conditional \U1D479\U00B2", "Individual", "Year", "Residual", "Individual", "Year", "Residual", "Years", "Individuals", "Observations"))# Fixed effect sizes (non-standardized)fixefTable <- stats_mod_max_score$tidy %>% dplyr::filter(effect == "fixed") %>% dplyr::select(term, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "fixed")# Fixed effect sizes (standardized)fixef_bw_Table <- stats_mod_max_score$partR2m$BW %>% as.data.frame() %>% mutate(stat = "fixed_bw") %>% rename(conf.low = CI_lower, conf.high = CI_upper)# Semi-partial R2 estimatesR2Table <- bind_rows(stats_mod_max_score$partR2m$R2, stats_mod_max_score$partR2c$R2[1,]) %>% dplyr::select(term, estimate, CI_lower, CI_upper) %>% as.data.frame() %>% mutate(stat = "partR2") %>% rename(conf.low = CI_lower, conf.high = CI_upper)# Random effects variancesranefTable <- stats_mod_max_score$tidy %>% dplyr::filter(effect == "ran_pars") %>% dplyr::select(group, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "rand") %>% rename(term = group) %>% mutate(estimate = estimate^2, conf.high = conf.high^2, conf.low = conf.low^2)# Adjusted repeatabilitiescoefRptTable <- stats_mod_max_score$rptR$R_boot %>% dplyr::select(-Fixed) %>% mutate(residual = 1 - rowSums(.)) %>% apply(., 2, function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% t() %>% as.data.frame() %>% rownames_to_column("term") %>% rename(estimate = V1, conf.low = `2.5%`, conf.high = `97.5%`) %>% mutate(stat = "RptR")# Store all parameters into a single table and clean it upallCoefs_mod <- bind_rows(fixef_bw_Table, R2Table, ranefTable, coefRptTable, sample_sizes) %>% bind_cols(., mod_comp_names) %>% mutate(coefString = ifelse(!is.na(conf.low), paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"), NA), effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)), rep("Partitioned \U1D479\U00B2", nrow(R2Table)), rep("Random effects \U1D70E\U00B2", nrow(ranefTable)), rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)), rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>% dplyr::select(effect, everything())# draw gt tablemod_max_score_table <- allCoefs_mod %>% dplyr::select(effect, comp_name, estimate, coefString) %>% gt(rowname_col = "row", groupname_col = "effect") %>% cols_label(comp_name = html("<i>Banded Dotterel rufous band score</i>"), estimate = "Mean estimate", coefString = "95% confidence interval") %>% fmt_number(columns = c(estimate), rows = 1:10, decimals = 2, use_seps = FALSE) %>% fmt_number(columns = c(estimate), rows = 11:13, decimals = 0, use_seps = FALSE) %>% sub_missing(columns = 1:4, missing_text = "") %>% cols_align(align = "left", columns = c(comp_name)) %>% tab_options(row_group.font.weight = "bold", row_group.background.color = brewer.pal(9,"Greys")[3], table.font.size = 12, data_row.padding = 3, row_group.padding = 4, summary_row.padding = 2, column_labels.font.size = 14, row_group.font.size = 12, table.width = pct(100))```Males tend to have a more full and immaculate rufous bands than females, however there is substantial overlap in their distributions.#### sex-specific breeding plumage distributions```{r, echo=FALSE, fig.width=6, fig.align='center'}ind_breeding_scores %>% group_by(rings_comb, sex) %>% summarise(max_breeding_score = max(max_breeding_score)) %>% ungroup() %>% group_by(sex, max_breeding_score) %>% summarise(n = n()) %>% bind_rows(., data.frame(sex = c("F", "M"), max_breeding_score = c(7, 3), n = c(0, 0))) %>% ungroup() %>% ggplot() + geom_bar(aes(x = max_breeding_score, y = n, fill = sex), alpha = 0.5, stat = "identity", position = "dodge", width = 0.5) + theme(text = element_text(family = "Verdana"), legend.title = element_blank(), axis.ticks = element_blank(), axis.text = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), axis.text.y = element_text(size = 10, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), strip.background = element_rect(fill = 'grey90', color = NA), strip.text = element_text(size = 13, colour = "grey40"), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_blank()) + theme(legend.position = c(0.15, 0.8)) + xlab("maximum individual rufous band score\nobserved during breeding season") + ylab("Frequency") + scale_colour_brewer(palette = "Dark2", direction = -1, name = "sex", labels = c("Female (N = 49)", "Male (N = 35)")) + scale_fill_brewer(palette = "Dark2", direction = -1, name = "sex", labels = c("Female (N = 49)", "Male (N = 35)"))```#### effect-size table for sex-specific breeding plumage model```{r, echo=FALSE}mod_max_score_table``````{r, include=FALSE}effects_mod_max_score <- as.data.frame(allEffects(mod_max_score)$sex) %>% mutate(x_mean = ifelse(sex == "F", as.numeric(factor(sex)) + 0.2, ifelse(sex == "M", as.numeric(factor(sex)) - 0.2, NA)), x_data = ifelse(sex == "F", as.numeric(factor(sex)), ifelse(sex == "M", as.numeric(factor(sex)), NA)))```#### sex-specific breeding plumage model predictions ```{r, echo=FALSE, fig.width=4, fig.align='center'}ggplot() + geom_errorbar(data = effects_mod_max_score, aes(x = x_mean, ymin = lower, ymax = upper), width = 0.05) + geom_point(data = effects_mod_max_score, aes(x = x_mean, y = fit, fill = sex), shape = 21, size = 5) + geom_jitter(data = ind_breeding_scores, aes(x = as.numeric(factor(sex)), y = max_breeding_score, fill = sex), width = 0.05, height = 0.1, shape = 21, size = 3, alpha = 0.75) + theme(legend.position = "none", text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), axis.text.y = element_text(size = 10, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), strip.background = element_rect(fill = 'grey90', color = NA), strip.text = element_text(size = 13, colour = "grey40"), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), axis.title.x = element_blank()) + scale_x_continuous(limits = c(0.8, 2.2), breaks = c(1.1, 1.9), labels = c("female", "male")) + ylab("maximum breeding plumage score observed\n(mean ± 95% CI)") + scale_fill_brewer(palette = "Set2") + scale_y_continuous(limits = c(3.8, 7.2))```### Age- and sex-specific variation in breeding plumage (of known-aged birds)data wrangle for sex-age model. Standardize max_breeding_score according to sex to make males and females to make them comparable across ages (i.e., does max_breeding_score change across age while accounting for underlying sex differences?)```{r}# extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual.ind_breeding_scores_age <- dat %>%mutate(score =as.numeric(score)) %>%mutate(breeding_season =ifelse(month(date) %in%c(8, 9, 10, 11, 12), 1, 0)) %>%filter(breeding_season ==1) %>%group_by(season, rings_comb, sex, age) %>%summarise(max_breeding_score =max(score)) %>%ungroup() %>%group_by(sex) %>%mutate(std_max_breeding_score = (max_breeding_score -mean(max_breeding_score, na.rm =TRUE)) /sd(max_breeding_score, na.rm =TRUE)) %>%filter(!is.na(age))```Assess sample sizes of each sex for individuals with known age```{r}ind_breeding_scores_age %>%group_by(sex) %>%summarise(n_distinct(rings_comb))```mixed model of standardized max breeding score (by sex) with individual as random intercept and age as fixed effect (Note: we tried to include season as an additional random intercept, but the low sample size of repeated measures across seasons caused convergence issues of the model)```{r}# linear mixed model for the difference in max score between sexes across agemod_max_score_age <-lmer(std_max_breeding_score ~ age + (1| rings_comb), data = ind_breeding_scores_age)``````{r, eval=FALSE, include=FALSE}# quick look at effectplot(allEffects(mod_max_score_age), type = "response")tbl_regression(mod_max_score_age, intercept = TRUE, label = list(age ~ "Age"))# Derive confidence intervals of effect sizes from parametric bootstrappingtidy_mod_max_score_age <- tidy(mod_max_score_age, conf.int = TRUE, conf.method = "boot", nsim = 1000)# run rptR to obtain repeatabilities of random effectsrpt_mod_max_score_age <- rpt(std_max_breeding_score ~ age + (1 | rings_comb), grname = c("rings_comb", "Fixed"), data = ind_breeding_scores_age, datatype = "Gaussian", nboot = 1000, npermut = 1000, ratio = TRUE, adjusted = TRUE, ncores = 4, parallel = TRUE)# run partR2 on each model to obtain marginal R2, parameter estimates, and beta# weightsR2m_mod_max_score_age <- partR2(mod_max_score_age, partvars = c("age"), R2_type = "marginal", nboot = 1000, CI = 0.95, max_level = 1)R2c_mod_max_score_age <- partR2(mod_max_score_age, partvars = c("age"), R2_type = "conditional", nboot = 1000, CI = 0.95, max_level = 1)stats_mod_max_score_age <- list(mod = mod_max_score_age, tidy = tidy_mod_max_score_age, rptR = rpt_mod_max_score_age, partR2m = R2m_mod_max_score_age, partR2c = R2c_mod_max_score_age, data = ind_breeding_scores_age)``````{r, include=FALSE}# saveRDS(stats_mod_max_score_age, file = "out/stats_mod_max_score_age.rds")stats_mod_max_score_age <- readRDS(file = "/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/out/stats_mod_max_score_age.rds")``````{r, include=FALSE}#### Table of effect sizes ----# Retrieve sample sizessample_sizes <- ind_breeding_scores_age %>% ungroup() %>% summarise(Year = n_distinct(season), Individual = n_distinct(rings_comb), Observations = nrow(.))sample_sizes <- as.data.frame(t(as.data.frame(sample_sizes))) %>% rownames_to_column("term") %>% rename(estimate = V1) %>% mutate(stat = "n")# clean model component namesmod_comp_names <- data.frame(comp_name = c("Age", "Total Marginal \U1D479\U00B2", "Age", "Total Conditional \U1D479\U00B2", "Individual", "Residual", "Individual", "Residual", "Years", "Individuals", "Observations"))# Fixed effect sizes (non-standardized)fixefTable <- stats_mod_max_score_age$tidy %>% dplyr::filter(effect == "fixed") %>% dplyr::select(term, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "fixed")# Fixed effect sizes (standardized)fixef_bw_Table <- stats_mod_max_score_age$partR2m$BW %>% as.data.frame() %>% mutate(stat = "fixed_bw") %>% rename(conf.low = CI_lower, conf.high = CI_upper)# Semi-partial R2 estimatesR2Table <- bind_rows(stats_mod_max_score_age$partR2m$R2, stats_mod_max_score_age$partR2c$R2[1,]) %>% dplyr::select(term, estimate, CI_lower, CI_upper) %>% as.data.frame() %>% mutate(stat = "partR2") %>% rename(conf.low = CI_lower, conf.high = CI_upper)# Random effects variancesranefTable <- stats_mod_max_score_age$tidy %>% dplyr::filter(effect == "ran_pars") %>% dplyr::select(group, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "rand") %>% rename(term = group) %>% mutate(estimate = estimate^2, conf.high = conf.high^2, conf.low = conf.low^2)# Adjusted repeatabilitiescoefRptTable <- stats_mod_max_score_age$rptR$R_boot %>% dplyr::select(-Fixed) %>% mutate(residual = 1 - rowSums(.)) %>% apply(., 2, function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% t() %>% as.data.frame() %>% rownames_to_column("term") %>% rename(estimate = V1, conf.low = `2.5%`, conf.high = `97.5%`) %>% mutate(stat = "RptR")# Store all parameters into a single table and clean it upallCoefs_mod <- bind_rows(fixef_bw_Table, R2Table, ranefTable, coefRptTable, sample_sizes) %>% bind_cols(., mod_comp_names) %>% mutate(coefString = ifelse(!is.na(conf.low), paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"), NA), effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)), rep("Partitioned \U1D479\U00B2", nrow(R2Table)), rep("Random effects \U1D70E\U00B2", nrow(ranefTable)), rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)), rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>% dplyr::select(effect, everything())# draw gt tablemod_max_score_age_table <- allCoefs_mod %>% dplyr::select(effect, comp_name, estimate, coefString) %>% gt(rowname_col = "row", groupname_col = "effect") %>% cols_label(comp_name = html("<i>Banded Dotterel rufous band score</i>"), estimate = "Mean estimate", coefString = "95% confidence interval") %>% fmt_number(columns = c(estimate), rows = 1:10, decimals = 2, use_seps = FALSE) %>% fmt_number(columns = c(estimate), rows = 9:11, decimals = 0, use_seps = FALSE) %>% sub_missing(columns = 1:4, missing_text = "") %>% cols_align(align = "left", columns = c(comp_name)) %>% tab_options(row_group.font.weight = "bold", row_group.background.color = brewer.pal(9,"Greys")[3], table.font.size = 12, data_row.padding = 3, row_group.padding = 4, summary_row.padding = 2, column_labels.font.size = 14, row_group.font.size = 12, table.width = pct(100))```#### effect-size table for sex- and age-specific variation in breeding plumageno effect of age on the maximum extent of an individual's rufous band score while controlling for underlying sex-specific variation```{r, echo=FALSE}mod_max_score_age_table```#### age-specific breeding plumage model predictions```{r, echo=FALSE}# extract predicted trendspr <- ggeffects::predict_response(mod_max_score_age)# join the back-transformed dates to model fitsmod_max_score_age_fits <- as.data.frame(pr)ggplot() + geom_ribbon(data = mod_max_score_age_fits, aes(x = age.x, ymin = age.conf.low, ymax = age.conf.high), fill = "grey70", alpha = 0.5) + geom_line(data = mod_max_score_age_fits, aes(x = age.x, y = age.predicted), color = "grey40") + geom_jitter(data = ind_breeding_scores_age, aes(x = age, y = std_max_breeding_score, fill = sex), width = 0.05, height = 0.1, shape = 21, size = 3, alpha = 0.75, alpha = 1) + theme(legend.position = c(0.15, 0.2), legend.title = element_blank(), text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), axis.text.y = element_text(size = 10, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), strip.background = element_rect(fill = 'grey90', color = NA), strip.text = element_text(size = 13, colour = "grey40"), panel.grid.major = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank()) + scale_x_continuous(limits = c(0.5, 6.5), breaks = c(1:6), labels = c(1:6)) + ylab("sex-standardized\nmaximum breeding plumage score observed") + xlab("age") + scale_fill_brewer(palette = "Dark2", direction = -1, name = "sex", labels = c("Female (N = 14)", "Male (N = 10)"))```### Migratory status- and sex-specific variation in breeding plumage (of known-status birds)data wrangle for sex-migratory status model```{r}# extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual.ind_breeding_scores_status <- dat %>%mutate(score =as.numeric(score)) %>%mutate(breeding_season =ifelse(month(date) %in%c(8, 9, 10, 11, 12), 1, 0)) %>%filter(breeding_season ==1) %>%group_by(season, rings_comb, sex, age, migratory_status) %>%summarise(max_breeding_score =max(score)) %>%ungroup() %>%group_by(sex) %>%mutate(std_max_breeding_score = (max_breeding_score -mean(max_breeding_score, na.rm =TRUE)) /sd(max_breeding_score, na.rm =TRUE)) %>%filter(migratory_status %in%c("M", "R"))```Assess sample sizes of each sex and migratory status (for birds with known status)```{r, echo=TRUE}# Assess sample sizes of each sex and season # (13 females (2 migrant, 11 resident), 8 males (3 migrant, 5 resident))ind_breeding_scores_status %>% group_by(sex, migratory_status) %>% summarise(n_distinct(rings_comb))```mixed model of standardized max breeding score (by sex) with individual as random intercept and migratory status as fixed effect (Note: we tried to include season as an additional random intercept, but the low sample size of repeated measures across seasons caused convergence issues of the model)```{r}# linear mixed model for the difference in max score between sexesmod_max_score_status <-lmer(std_max_breeding_score ~ migratory_status + (1| rings_comb), data = ind_breeding_scores_status)``````{r, eval=FALSE, include=FALSE}# # detected singularity issue in a model with (1 | season), summary showed # that variance component of this random effect was zero. # # Model with one random effectmod1 <- lmer(std_max_breeding_score ~ migratory_status + season + (1 | rings_comb), data = ind_breeding_scores_status)mod2 <- lmer(std_max_breeding_score ~ migratory_status + (1 | season), data = ind_breeding_scores_status)mod3 <- lmer(std_max_breeding_score ~ migratory_status + (1 | rings_comb), data = ind_breeding_scores_status)# Compare with the full modelmod_full <- lmer(std_max_breeding_score ~ migratory_status + (1 | rings_comb) + (1 | season), data = ind_breeding_scores_status)# Check the summariessummary(mod1)summary(mod2)summary(mod3)summary(mod_full)# Compare models using AICAIC(mod1, mod2, mod3, mod_full)# conclude that the term (1 | season) is not needed for statistical reasons, but also since we are not interested in population-level annual variation in the max-breeding plumage scores we can drop this``````{r, eval=FALSE, include=FALSE}# quick look at effectplot(allEffects(mod_max_score_status), type = "response")tbl_regression(mod_max_score_status, intercept = TRUE, label = list(migratory_status ~ "migratory status"))# Derive confidence intervals of effect sizes from parametric bootstrappingtidy_mod_max_score_status <- tidy(mod_max_score_status, conf.int = TRUE, conf.method = "boot", nsim = 1000)# run rptR to obtain repeatabilities of random effectsrpt_mod_max_score_status <- rpt(std_max_breeding_score ~ migratory_status + (1 | rings_comb), grname = c("rings_comb", "Fixed"), data = ind_breeding_scores_status, datatype = "Gaussian", nboot = 1000, npermut = 1000, ratio = TRUE, adjusted = TRUE, ncores = 4, parallel = TRUE)# run partR2 on each model to obtain marginal R2, parameter estimates, and beta# weightsR2m_mod_max_score_status <- partR2(mod_max_score_status, partvars = c("migratory_status"), R2_type = "marginal", nboot = 1000, CI = 0.95, max_level = 1)R2c_mod_max_score_status <- partR2(mod_max_score_status, partvars = c("migratory_status"), R2_type = "conditional", nboot = 1000, CI = 0.95, max_level = 1)stats_mod_max_score_status <- list(mod = mod_max_score_status, tidy = tidy_mod_max_score_status, rptR = rpt_mod_max_score_status, partR2m = R2m_mod_max_score_status, partR2c = R2c_mod_max_score_status, data = ind_breeding_scores_status)``````{r, include=FALSE}# saveRDS(stats_mod_max_score_status, file = "out/stats_mod_max_score_status.rds")stats_mod_max_score_status <- readRDS(file = "/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/out/stats_mod_max_score_status.rds")``````{r, echo=FALSE}#### Table of effect sizes ----# Retrieve sample sizessample_sizes <- ind_breeding_scores_status %>% ungroup() %>% summarise(Year = n_distinct(season), Individual = n_distinct(rings_comb), Observations = nrow(.))sample_sizes <- as.data.frame(t(as.data.frame(sample_sizes))) %>% rownames_to_column("term") %>% rename(estimate = V1) %>% mutate(stat = "n")# clean model component namesmod_comp_names <- data.frame(comp_name = c("Resident (Migratory Status)", "Total Marginal \U1D479\U00B2", "Migratory Status", "Total Conditional \U1D479\U00B2", "Individual", "Residual", "Individual", "Residual", "Years", "Individuals", "Observations"))# Fixed effect sizes (non-standardized)fixefTable <- stats_mod_max_score_status$tidy %>% dplyr::filter(effect == "fixed") %>% dplyr::select(term, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "fixed")# Fixed effect sizes (standardized)fixef_bw_Table <- stats_mod_max_score_status$partR2m$BW %>% as.data.frame() %>% mutate(stat = "fixed_bw") %>% rename(conf.low = CI_lower, conf.high = CI_upper)# Semi-partial R2 estimatesR2Table <- bind_rows(stats_mod_max_score_status$partR2m$R2, stats_mod_max_score_status$partR2c$R2[1,]) %>% dplyr::select(term, estimate, CI_lower, CI_upper) %>% as.data.frame() %>% mutate(stat = "partR2") %>% rename(conf.low = CI_lower, conf.high = CI_upper)# Random effects variancesranefTable <- stats_mod_max_score_status$tidy %>% dplyr::filter(effect == "ran_pars") %>% dplyr::select(group, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "rand") %>% rename(term = group) %>% mutate(estimate = estimate^2, conf.high = conf.high^2, conf.low = conf.low^2)# Adjusted repeatabilitiescoefRptTable <- stats_mod_max_score_status$rptR$R_boot %>% dplyr::select(-Fixed) %>% mutate(residual = 1 - rowSums(.)) %>% apply(., 2, function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% t() %>% as.data.frame() %>% rownames_to_column("term") %>% rename(estimate = V1, conf.low = `2.5%`, conf.high = `97.5%`) %>% mutate(stat = "RptR")# Store all parameters into a single table and clean it upallCoefs_mod <- bind_rows(fixef_bw_Table, R2Table, ranefTable, coefRptTable, sample_sizes) %>% bind_cols(., mod_comp_names) %>% mutate(coefString = ifelse(!is.na(conf.low), paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"), NA), effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)), rep("Partitioned \U1D479\U00B2", nrow(R2Table)), rep("Random effects \U1D70E\U00B2", nrow(ranefTable)), rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)), rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>% dplyr::select(effect, everything())mod_max_score_table_status <- allCoefs_mod %>% dplyr::select(effect, comp_name, estimate, coefString) %>% gt(rowname_col = "row", groupname_col = "effect") %>% cols_label(comp_name = html("<i>Banded Dotterel rufous band score</i>"), estimate = "Mean estimate", coefString = "95% confidence interval") %>% fmt_number(columns = c(estimate), rows = 1:8, decimals = 2, use_seps = FALSE) %>% fmt_number(columns = c(estimate), rows = 9:11, decimals = 0, use_seps = FALSE) %>% sub_missing(columns = 1:4, missing_text = "") %>% cols_align(align = "left", columns = c(comp_name)) %>% tab_options(row_group.font.weight = "bold", row_group.background.color = brewer.pal(9,"Greys")[3], table.font.size = 12, data_row.padding = 3, row_group.padding = 4, summary_row.padding = 2, column_labels.font.size = 14, row_group.font.size = 12, table.width = pct(80))```#### effect-size table and plot for sex- and migratory-status breeding plumageno effect of migratory status on an individual's rufous band score while controlling for underlying sex-specific variation```{r, echo=FALSE}mod_max_score_table_status```#### migratory status-specific breeding plumage model predictions ```{r, echo=FALSE, fig.width=6, fig.align='center'}effects_mod_max_score_status2 <- as.data.frame(allEffects(mod_max_score_status)$migratory_status) %>% mutate(x_mean = ifelse(migratory_status == "M", as.numeric(factor(migratory_status)) + 0.2, ifelse(migratory_status == "R", as.numeric(factor(migratory_status)) - 0.2, NA)), x_data = ifelse(migratory_status == "M", as.numeric(factor(migratory_status)), ifelse(migratory_status == "R", as.numeric(factor(migratory_status)), NA)))ggplot() + geom_errorbar(data = effects_mod_max_score_status2, aes(x = x_mean, ymin = lower, ymax = upper), width = 0.05) + geom_point(data = effects_mod_max_score_status2, aes(x = x_mean, y = fit, fill = migratory_status), shape = 21, size = 5) + geom_jitter(data = ind_breeding_scores_status, aes(x = as.numeric(factor(migratory_status)), y = std_max_breeding_score, fill = migratory_status), width = 0.05, height = 0.1, shape = 21, size = 3, alpha = 0.75) + theme(legend.position = "none", text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), axis.title.x = element_blank()) + scale_x_continuous(limits = c(0.8, 2.2), breaks = c(1.1, 1.9), labels = c("migrant", "resident")) + ylab("sex-standardized\nmaximum breeding plumage score observed") + scale_fill_brewer(palette = "Set1")```## Sources of individual variation in the timing of pre-basic moult### Within-individual moult dynamics across sexeswrangle data by calculating the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart sex-specific variation in the timing of pre-basic moult```{r, include=FALSE}ind_prop_molt_scores <- dat %>% mutate(score = as.numeric(score)) %>% left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% filter(!is.na(max_breeding_score)) %>% filter(date_J_of_max_score <= date_J) %>% group_by(sex) %>% mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>% # focus on key months of moult (Nov, Dec, Jan, Feb, March) filter(month(date) %in% c(11, 12, 1, 2, 3)) %>% # subtract 1 to make minimum score 0 mutate(prop_molt_score_max = (score-1)/(max_breeding_score-1), prop_molt_score_mode = (score-1)/(mode_breeding_score-1), prop_molt_score_max_mean = (score-1)/(max_breeding_score_mean-1)) %>% select(rings_comb, season, date_J, sex, migratory_status, age, score, max_breeding_score, mode_breeding_score, max_breeding_score_mean, std_max_breeding_score, prop_molt_score_max, prop_molt_score_mode, prop_molt_score_max_mean) %>% arrange(rings_comb, season, date_J) %>% distinct() %>% filter(!is.na(prop_molt_score_max)) %>% ungroup() %>% mutate(prop_molt_score_max_mean = ifelse(prop_molt_score_max_mean > 1, 1, prop_molt_score_max_mean), prop_molt_score_mode = ifelse(prop_molt_score_mode > 1, 1, prop_molt_score_mode)) %>% mutate(ring_season = paste(rings_comb, season, sep = "_"))```Assess sample sizes of each sex```{r, echo=FALSE}ind_prop_molt_scores %>% group_by(sex) %>% summarise(n_distinct(rings_comb))```mixed effects binomial model assessing sex-differences in the timing of pre-basic moult while accounting for bird- and year-specific random variation in the timing of pre-basic moult (i.e., random slope model of date for individual and year).```{r}mod_moult_sex_bird_RIS <- lme4::glmer(prop_molt_score_max ~ date_J + sex + (date_J | ring_season),data = ind_prop_molt_scores, family = binomial)``````{r, echo=FALSE}tbl_regression(mod_moult_sex_bird_RIS, intercept = TRUE, label = list(date_J ~ "Date", sex ~ "Sex"))``````{r, echo=FALSE}RIS_bird_sex <- augment(mod_moult_sex_bird_RIS) %>% dplyr::select(prop_molt_score_max, date_J, ring_season, .fitted, sex)dates_for_plot <- data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), date_J = c(1:365))df_bird_sex <- RIS_bird_sex %>% dplyr::group_by(sex, ring_season, date_J) %>% dplyr::summarise(RIS = mean(.fitted)) %>% mutate(RIS_trans = exp(RIS)/(1 + exp(RIS))) %>% left_join(., dates_for_plot, by = "date_J") %>% mutate(rings_comb = as.factor(ring_season)) %>% filter(date > as.Date("2021-11-01"))# extract predicted trendspr <- ggeffects::predict_response(mod_moult_sex_bird_RIS, c("date_J [30:293]", "sex"))# join the back-transformed dates to model fitsmod_moult_sex_bird_RIS_fits <- as.data.frame(pr) %>% rename(date_J = x, sex = group) %>% left_join(., dates_for_plot, by = "date_J") %>% filter(date > as.Date("2021-11-01"))plot_rs_bird_sex <- ggplot() + geom_line(data = df_bird_sex, aes(x = date, y = RIS_trans, group = ring_season), alpha = 0.25, colour = "grey50") + # geom_line(data = df_bird_season %>% filter(rings_comb == "RRGO"), #"RBRL"), # aes(x = date, y = RIS_trans, # group = bird_season), # alpha = 1, colour = "black") + # geom_jitter(data = mid_point_moult_df, # aes(x = date, y = 0.75), # height = 0.025, color = "red") + # geom_jitter(data = moult_commence_ring_season, # aes(x = date, y = 0.25), # height = 0.025, color = "blue") + # geom_point(data = mid_point_moult_df %>% filter(rings_comb == "RRGO"), # aes(x = date, y = 0.55)) + # geom_point(data = moult_commence_ring_season %>% filter(rings_comb == "RRGO"), # aes(x = date, y = 0.45)) + theme(legend.position = "none", legend.justification = c(1, 0), legend.title = element_blank(), text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), axis.text.y = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + scale_y_continuous(limits = c(-0.05, 1.05)) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month") + ggtitle("random slopes of\npre-basic moult timing")# plot the modelplot_rs_bird_fs_sex <- ggplot() + geom_line(data = mod_moult_sex_bird_RIS_fits, aes(x = date, y = predicted, color = sex)) + geom_ribbon(data = mod_moult_sex_bird_RIS_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = sex), lwd = 1, alpha = 0.25) + geom_jitter(data = ind_prop_molt_scores %>% left_join(., dates_for_plot, by = "date_J"), aes(x = date, y = prop_molt_score_max, color = sex), alpha = 0.25, shape = 20, width = 0, height = 0.025) + theme(legend.position = c(0.55, 0.1), # legend.background = element_rect(fill = 'transparent', color = NA), legend.justification = c(1, 0), legend.title = element_blank(), text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1, colour = "grey40"), axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + scale_y_continuous(limits = c(-0.05, 1.05)) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month") + scale_colour_brewer(palette = "Dark2", direction = -1, name = "Sex", labels = c("Female (N = 49)", "Male (N = 32)")) + scale_fill_brewer(palette = "Dark2", direction = -1, name = "Sex", labels = c("Female (N = 49)", "Male (N = 32)")) + ggtitle("sex-specific slopes\nof pre-basic moult timing")plot_rs_bird_sex + plot_rs_bird_fs_sex ``````{r, include=FALSE, eval=FALSE}# # quick look at effect# plot(allEffects(mod_moult_sex_bird_season_RIS), type = "response")# tbl_regression(mod_moult_sex_bird_season_RIS)# # # not much variation in slopes across seasons, so simplify the model to only include the (date_J | rings_comb) term.# mod_moult_sex_bird_RIS <- # lme4::glmer(prop_molt_score ~ # date_J + sex + # (date_J | rings_comb),# data = ind_prop_molt_scores, # family = binomial)# # tab_model(mod_moult_sex_bird_RIS)# # # Load the package# library(performance)# # # Calculate R-squared# r2(mod_moult_sex_bird_RIS)# # # Load the package# library(MuMIn)# # # Calculate R-squared# r.squaredGLMM(mod_moult_sex_bird_RIS)# # # Derive confidence intervals of effect sizes from parametric bootstrapping# tidy_mod_moult_sex_bird_RIS <-# tidy(mod_moult_sex_bird_RIS, conf.int = TRUE, conf.method = "boot", nsim = 1000)# # saveRDS(tidy_mod_moult_sex_bird_RIS, here("out/tidy_mod_moult_sex_bird_RIS.rds"))``````{r, include=FALSE}# derive the mid-point date of each individual's moult# wrangle the dates of plumage retention and moult for subsequent analysis# Step 1: Calculate the max_score for each individual and seasonmax_scores <- ind_prop_molt_scores %>% group_by(rings_comb, season) %>% summarise(max_score = max(prop_molt_score_max, na.rm = TRUE), .groups = "drop")# Step 2: Merge the max_scores back into the original datamid_point_dates <- ind_prop_molt_scores %>% left_join(max_scores, by = c("rings_comb", "season")) %>% # Add the max_score group_by(rings_comb, season) %>% arrange(date_J) %>% # Ensure data is sorted by date_J # Step 3: Find the latest date_J for the maximum prop_molt_score filter(prop_molt_score_max == max_score) %>% slice_max(order_by = date_J, n = 1, with_ties = FALSE) %>% rename(latest_max_date_J = date_J) %>% # Step 4: Find the earliest date_J where prop_molt_score < max_score right_join( ind_prop_molt_scores %>% left_join(max_scores, by = c("rings_comb", "season")) %>% # Add the max_score here too group_by(rings_comb, season) %>% arrange(date_J) %>% filter(prop_molt_score_max < max_score) %>% slice_min(order_by = date_J, n = 1, with_ties = FALSE) %>% rename(earliest_decline_date_J = date_J), by = c("rings_comb", "season") ) %>% # Step 5: Compute the mid-point date_J mutate(mid_point_date_J = round((latest_max_date_J + earliest_decline_date_J) / 2)) %>% # Select relevant columns select(rings_comb, season, latest_max_date_J, earliest_decline_date_J, mid_point_date_J) %>% distinct()dates_for_plot <- data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), date_J = c(1:365))moult_commence_ring_season <- left_join(mid_point_dates, dates_for_plot %>% rename(mid_point_date_J = date_J), by = "mid_point_date_J") %>% mutate(season = as.factor(season))```### Within-individual moult dynamics across migratory status```{r, include=FALSE}# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status.ind_prop_molt_scores_status <- ind_prop_molt_scores %>% filter(migratory_status %in% c("R", "M"))``````{r, echo=FALSE}# Assess sample sizes of each migratory statusind_prop_molt_scores_status %>% group_by(migratory_status) %>% summarise(n_distinct(rings_comb))``````{r}# mixed effects binomial model comparing migratory status and date effect on the changes in moult scoresmod_moult_status_bird_RIS <- lme4::glmer(prop_molt_score_max ~ date_J + migratory_status + (date_J|ring_season),data = ind_prop_molt_scores_status, family = binomial)``````{r, echo=FALSE}tbl_regression(mod_moult_status_bird_RIS, intercept = TRUE, label = list(date_J ~ "Date", migratory_status ~ "Migratory status"))``````{r, echo=FALSE}RIS_bird_status <- augment(mod_moult_status_bird_RIS) %>% dplyr::select(prop_molt_score_max, date_J, ring_season, .fitted, migratory_status)dates_for_plot <- data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), date_J = c(1:365))df_bird_status <- RIS_bird_status %>% dplyr::group_by(migratory_status, ring_season, date_J) %>% dplyr::summarise(RIS = mean(.fitted)) %>% mutate(RIS_trans = exp(RIS)/(1 + exp(RIS))) %>% left_join(., dates_for_plot, by = "date_J") %>% mutate(rings_comb = as.factor(ring_season)) %>% filter(date > as.Date("2021-11-01"))# extract predicted trendspr <- ggeffects::predict_response(mod_moult_status_bird_RIS, c("date_J [30:293]", "migratory_status"))# join the back-transformed dates to model fitsmod_moult_status_bird_RIS_fits <- as.data.frame(pr) %>% rename(date_J = x, migratory_status = group) %>% left_join(., dates_for_plot, by = "date_J") %>% filter(date > as.Date("2021-11-01"))plot_rs_bird_status <- ggplot() + geom_line(data = df_bird_status, aes(x = date, y = RIS_trans, group = ring_season), alpha = 0.25, colour = "grey50") + theme(legend.position = "none", legend.justification = c(1, 0), legend.title = element_blank(), text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), axis.text.y = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + scale_y_continuous(limits = c(-0.05, 1.05)) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month") + ggtitle("bird-specific slopes of\npre-basic moult timing (random)")# plot the modelplot_rs_bird_status_fs <- ggplot() + geom_line(data = mod_moult_status_bird_RIS_fits, aes(x = date, y = predicted, color = migratory_status)) + geom_ribbon(data = mod_moult_status_bird_RIS_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = migratory_status), lwd = 1, alpha = 0.25) + geom_jitter(data = ind_prop_molt_scores_status %>% left_join(., dates_for_plot, by = "date_J"), aes(x = date, y = prop_molt_score_max, color = migratory_status), alpha = 0.25, shape = 20, width = 0, height = 0.025) + theme(legend.position = c(0.55, 0.1), # legend.background = element_rect(fill = 'transparent', color = NA), legend.justification = c(1, 0), legend.title = element_blank(), text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1, colour = "grey40"), axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + scale_y_continuous(limits = c(-0.05, 1.05)) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month") + scale_colour_brewer(palette = "Set1", labels = c("Migrant (N = 5)", "Resident (N = 16)")) + scale_fill_brewer(palette = "Set1", labels = c("Migrant (N = 5)", "Resident (N = 16)")) + ggtitle("status-specific slopes\nof pre-basic moult timing (fixed)")plot_rs_bird_status + plot_rs_bird_status_fs```### Within-individual moult dynamics across ages```{r, include=FALSE}# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation between age groups.ind_prop_molt_scores_age <- ind_prop_molt_scores %>% filter(!is.na(age))``````{r, echo=FALSE}# Assess sample sizes of each sexind_prop_molt_scores_age %>% group_by(age) %>% summarise(n_distinct(rings_comb))``````{r}# mixed effects binomial model comparing migratory status and date effect on the changes in moult scoresmod_moult_age_bird_RIS <- lme4::glmer(prop_molt_score_max ~ date_J + age + (date_J | ring_season),data = ind_prop_molt_scores_age, family = binomial)``````{r, include=FALSE}# ind_prop_molt_scores_age %>% # mutate(id_year = paste(rings_comb, season, sep = "_")) %>% # group_by(id_year) %>% # summarise(n_ = n()) %>% View()``````{r,echo=FALSE}tbl_regression(mod_moult_age_bird_RIS, intercept = TRUE, label = list(date_J ~ "Date", age ~ "Age"))``````{r, echo=FALSE}RIS_bird_age <- augment(mod_moult_age_bird_RIS) %>% dplyr::select(prop_molt_score_max, date_J, ring_season, .fitted, age)dates_for_plot <- data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), date_J = c(1:365))df_bird_age <- RIS_bird_age %>% dplyr::group_by(age, ring_season, date_J) %>% dplyr::summarise(RIS = mean(.fitted)) %>% mutate(RIS_trans = exp(RIS)/(1 + exp(RIS))) %>% left_join(., dates_for_plot, by = "date_J") %>% mutate(id_year = as.factor(ring_season)) %>% filter(date > as.Date("2021-11-01"))# extract predicted trendspr <- ggeffects::predict_response(mod_moult_age_bird_RIS, c("date_J [30:293]", "age"))# join the back-transformed dates to model fitsmod_moult_age_bird_RIS_fits <- as.data.frame(pr) %>% rename(date_J = x, age = group) %>% left_join(., dates_for_plot, by = "date_J") %>% filter(date > as.Date("2021-11-01"))plot_rs_bird_age <- ggplot() + geom_line(data = df_bird_age, aes(x = date, y = RIS_trans, group = ring_season), alpha = 0.25, colour = "grey50") + theme(legend.position = "none", legend.justification = c(1, 0), legend.title = element_blank(), text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), axis.text.y = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + scale_y_continuous(limits = c(-0.05, 1.05)) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month") + ggtitle("bird-specific slopes of\npre-basic moult timing (random)")# plot the modelplot_rs_bird_age_fs <- ggplot() + geom_line(data = mod_moult_age_bird_RIS_fits %>% mutate(age = as.factor(age)), aes(x = date, y = predicted, color = age)) + geom_ribbon(data = mod_moult_age_bird_RIS_fits %>% mutate(age = as.factor(age)), aes(x = date, ymax = conf.high, ymin = conf.low, fill = age), lwd = 1, alpha = 0.25) + geom_jitter(data = ind_prop_molt_scores_age %>% left_join(., dates_for_plot, by = "date_J") %>% mutate(age = as.factor(age)), aes(x = date, y = prop_molt_score_max, color = age), alpha = 1, shape = 20, width = 0, height = 0.025) + theme(legend.position = c(0.25, 0.1), # legend.background = element_rect(fill = 'transparent', color = NA), legend.justification = c(1, 0), legend.title = element_blank(), text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1, colour = "grey40"), axis.text.y = element_blank(), axis.title.y = element_blank(), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + scale_y_continuous(limits = c(-0.05, 1.05)) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month") + scale_colour_brewer(palette = "Spectral", direction = -1, name = "Age", labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", "4 (N = 13", "5 (N = 4)", "6 (N = 2)")) + scale_fill_brewer(palette = "Spectral", direction = -1, name = "Age", labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", "4 (N = 13", "5 (N = 4)", "6 (N = 2)")) + ggtitle("age-specific slopes\nof pre-basic moult timing (fixed)")plot_rs_bird_age + plot_rs_bird_age_fs```inspect the random slopes against the raw data from photos```{r, echo=FALSE, fig.height=20}ggplot() + geom_line(data = df_bird_age %>% mutate(age = as.factor(age)), aes(x = date, y = RIS_trans),#, color = age), alpha = 1) + geom_point(data = ind_prop_molt_scores_age %>% left_join(., dates_for_plot, by = "date_J") %>% mutate(age = as.factor(age), id_year = paste(rings_comb, season, sep = "_")), aes(x = date, y = prop_molt_score_max),# color = age), alpha = 1, shape = 20) + theme(legend.position = "none", legend.justification = c(1, 0), legend.title = element_blank(), text = element_text(family = "Verdana"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), axis.text.y = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + scale_y_continuous(limits = c(-0.05, 1.05)) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month") + facet_wrap("id_year", ncol = 3) #+ # scale_colour_brewer(palette = "Spectral", direction = -1)```## Breeding status analysis```{r, include=FALSE}# explore datasets provided by Ted## 2021 season# import and consolidate key columnsbreeding_data_2021 <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_2021_Feb2022_upd_Nov22_Locations.xlsx", sheet = "Kaikoura_NestVisit2021_1", col_types = "text") %>% mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>% mutate(date = as.Date(visit_date_nz)) %>% select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>% rename(nest_fate = `Nest Status`, eggs = `Egg count`, hatched = `Number Hatched`, fledged = `Number Fledged`, chicks = `Number of chicks seen`)# Define a function to extract the desired patterns and ensure enough columnsextract_patterns <- function(text) { # Replace vertical bar '|' with an empty string text <- gsub("\\|", "", text) # Extract all 4 or 5-character long texts that start with R or r matches_r <- str_extract_all(text, "\\b[Rr]\\w{3,4}\\b")[[1]] # Extract UB, UN, UBF, UBM as standalone patterns matches_ub_un <- str_extract_all(text, "\\b[Uu][BbNnFfMm]\\b")[[1]] # Combine matches matches <- c(matches_r, matches_ub_un) # Ensure all matches are capitalized matches <- toupper(matches) # Ensure there are exactly two columns, filling with NA if necessary result <- c(matches, rep(NA, 2 - length(matches))) # Return the result as a named vector return(setNames(result, c("parent1", "parent2")))}# Apply the function to the text column and store the results in new columnsextracted_parents <- t(sapply(breeding_data_2021$NestID, extract_patterns))breeding_data_2021 <- bind_cols(breeding_data_2021, extracted_parents) %>% select(-parent2) %>% rename(parent = parent1) %>% bind_rows(bind_cols(breeding_data_2021, extracted_parents) %>% select(-parent1) %>% rename(parent = parent2)) %>% filter(!is.na(parent)) %>% filter(parent %in% (dat %>% filter(season == 2021) %>% pull(rings_comb) %>% unique())) %>% filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>% group_by(parent) %>% mutate( first_date_breeding = date[which.min(date)], last_date_breeding = date[which.max(date)]) %>% # arrange(desc(date)) %>% # slice(1) %>% ungroup() %>% # arrange(desc(date)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>% # rename(last_date_breeding = date) %>% select(parent, first_date_breeding, last_date_breeding, nest_id, OBJECTID)# import and consolidate key columnsbreeding_data_2021_ <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_2021_Feb2022_upd_Nov22_Locations.xlsx", sheet = "Kaikoura_DotterelNest2021_0", col_types = "text") %>% mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"), fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"), found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"), fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"), found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>% mutate(hatch_date = as.Date(hatch_date_nz), fail_date = as.Date(fail_date_nz), found_date = as.Date(found_date_nz)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>% select(OBJECTID, nest_id, found_date, hatch_date, fail_date)breeding_2021_final <- left_join(breeding_data_2021, breeding_data_2021_, by = "nest_id") %>% mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) | ((!is.na(hatch_date) & found_date < hatch_date)) | ((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) | ((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>% group_by(parent) %>% mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% mutate(first_date_breeding_ = min(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% mutate(days_diff = last_date_breeding_ - last_date_breeding) %>% mutate(days_diff = first_date_breeding_ - first_date_breeding) %>% arrange(desc(days_diff)) %>% mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding)), first_date_breeding_final = as.Date(ifelse(days_diff > 0, first_date_breeding + ceiling(days_diff/2), first_date_breeding))) %>% # arrange(desc(last_date_breeding_final)) %>% select(parent, first_date_breeding_final, last_date_breeding_final) %>% # specify the season as the first calender year mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final))) %>% distinct()############# 2022 season# import and consolidate key columnsbreeding_data_2022 <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Jan23_upd.xlsx", sheet = "Kaikoura_NestVisit_1", col_types = "text") %>% mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>% mutate(date = as.Date(visit_date_nz)) %>% select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>% rename(nest_fate = `Nest Status`, eggs = `Egg count`, hatched = `Number Hatched`, fledged = `Number Fledged`, chicks = `Number of chicks seen`)# Apply the function to the text column and store the results in new columnsextracted_parents2 <- t(sapply(breeding_data_2022$NestID, extract_patterns))breeding_data_2022 <- bind_cols(breeding_data_2022, extracted_parents2) %>% select(-parent2) %>% rename(parent = parent1) %>% bind_rows(bind_cols(breeding_data_2022, extracted_parents2) %>% select(-parent1) %>% rename(parent = parent2)) %>% filter(!is.na(parent)) %>% filter(parent %in% (dat %>% filter(season == 2022) %>% pull(rings_comb) %>% unique())) %>% filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>% group_by(parent) %>% mutate( first_date_breeding = date[which.min(date)], last_date_breeding = date[which.max(date)]) %>% # arrange(desc(date)) %>% # slice(1) %>% ungroup() %>% # arrange(desc(date)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>% # rename(last_date_breeding = date) %>% select(parent, first_date_breeding, last_date_breeding, nest_id, OBJECTID)# import and consolidate key columnsbreeding_data_2022_ <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Jan23_upd.xlsx", sheet = "Kaikoura_DotterelNest_0", col_types = "text") %>% mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"), fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"), found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"), fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"), found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>% mutate(hatch_date = as.Date(hatch_date_nz), fail_date = as.Date(fail_date_nz), found_date = as.Date(found_date_nz)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>% select(OBJECTID, nest_id, found_date, hatch_date, fail_date)breeding_2022_final <- left_join(breeding_data_2022, breeding_data_2022_, by = "nest_id") %>% mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) | ((!is.na(hatch_date) & found_date < hatch_date)) | ((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) | ((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>% group_by(parent) %>% mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% mutate(first_date_breeding_ = min(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% mutate(days_diff = last_date_breeding_ - last_date_breeding) %>% mutate(days_diff = first_date_breeding_ - first_date_breeding) %>% arrange(desc(days_diff)) %>% mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding)), first_date_breeding_final = as.Date(ifelse(days_diff > 0, first_date_breeding + ceiling(days_diff/2), first_date_breeding))) %>% # arrange(desc(last_date_breeding_final)) %>% select(parent, first_date_breeding_final, last_date_breeding_final) %>% # specify the season as the first calender year mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final))) %>% distinct()############# 2023 season# import and consolidate key columnsbreeding_data_2023 <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Feb2024_upd.xlsx", sheet = "Kaikoura_NestVisit_1", col_types = "text") %>% mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>% mutate(date = as.Date(visit_date_nz)) %>% select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>% rename(nest_fate = `Nest Status`, eggs = `Egg count`, hatched = `Number Hatched`, fledged = `Number Fledged`, chicks = `Number of chicks seen`)# Apply the function to the text column and store the results in new columnsextracted_parents3 <- t(sapply(breeding_data_2023$NestID, extract_patterns))breeding_data_2023 <- bind_cols(breeding_data_2023, extracted_parents3) %>% select(-parent2) %>% rename(parent = parent1) %>% bind_rows(bind_cols(breeding_data_2023, extracted_parents3) %>% select(-parent1) %>% rename(parent = parent2)) %>% filter(!is.na(parent)) %>% filter(parent %in% (dat %>% filter(season == 2023) %>% pull(rings_comb) %>% unique())) %>% filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>% group_by(parent) %>% mutate( first_date_breeding = date[which.min(date)], last_date_breeding = date[which.max(date)]) %>% # arrange(desc(date)) %>% # slice(1) %>% ungroup() %>% # arrange(desc(date)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>% # rename(last_date_breeding = date) %>% select(parent, first_date_breeding, last_date_breeding, nest_id, OBJECTID)# import and consolidate key columnsbreeding_data_2023_ <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Feb2024_upd.xlsx", sheet = "Kaikoura_DotterelNest_0", col_types = "text") %>% mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"), fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"), found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"), fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"), found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>% mutate(hatch_date = as.Date(hatch_date_nz), fail_date = as.Date(fail_date_nz), found_date = as.Date(found_date_nz)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>% select(OBJECTID, nest_id, found_date, hatch_date, fail_date)breeding_2023_final <- left_join(breeding_data_2023, breeding_data_2023_, by = "nest_id") %>% mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) | ((!is.na(hatch_date) & found_date < hatch_date)) | ((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) | ((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>% group_by(parent) %>% mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% mutate(first_date_breeding_ = min(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% mutate(days_diff = last_date_breeding_ - last_date_breeding) %>% mutate(days_diff = first_date_breeding_ - first_date_breeding) %>% arrange(desc(days_diff)) %>% mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding)), first_date_breeding_final = as.Date(ifelse(days_diff > 0, first_date_breeding + ceiling(days_diff/2), first_date_breeding))) %>% # arrange(desc(last_date_breeding_final)) %>% select(parent, first_date_breeding_final, last_date_breeding_final) %>% # specify the season as the first calender year mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final))) %>% distinct()# combine breeding data from the 3 seasons into one tablebreeding_data_all <- bind_rows(breeding_2021_final,breeding_2022_final,breeding_2023_final) %>% rename(rings_comb = parent) %>% # change to Julian date shifted for the Southern Hemisphere (1 = July 1) mutate(last_breeding_date_J = as.numeric(format(last_date_breeding_final + 181, "%j")), first_breeding_date_J = as.numeric(format(first_date_breeding_final + 181, "%j"))) %>% ungroup() %>% # standardize the last breeding date by year mutate(last_breeding_date_J_s = scale_by(last_breeding_date_J ~ season, ., scale = 0), first_breeding_date_J_s = scale_by(first_breeding_date_J ~ season, ., scale = 0)) %>% # make the scaled date variable numeric class mutate(last_breeding_date_J_snum = as.numeric(last_breeding_date_J_s), first_breeding_date_J_snum = as.numeric(first_breeding_date_J_s))```### Relationships between breeding schedule and plumage```{r, include=FALSE}ind_prop_molt_scores_last_breeding_date <- ind_prop_molt_scores %>% full_join(breeding_data_all, by = c("rings_comb", "season")) %>% mutate(date_diff = abs(last_breeding_date_J - date_J)) %>% arrange(rings_comb, season, date_diff) %>% group_by(rings_comb, season) %>% slice(1) %>% filter(date_diff < 28) %>% # duration of incubation (days) dplyr::select(rings_comb, season, sex, last_breeding_date_J, last_breeding_date_J_snum, prop_molt_score_max, max_breeding_score, age, migratory_status) %>% distinct() %>% group_by(sex) %>% mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>% ungroup() %>% left_join(., moult_commence_ring_season %>% mutate(season = as.numeric(as.character(season))), by = c("rings_comb", "season"))ind_prop_molt_scores_first_breeding_date <- ind_prop_molt_scores %>% full_join(breeding_data_all, by = c("rings_comb", "season")) %>% mutate(date_diff = abs(first_breeding_date_J - date_J)) %>% arrange(rings_comb, season, date_diff) %>% group_by(rings_comb, season) %>% slice(1) %>% dplyr::select(rings_comb, season, sex, first_breeding_date_J, first_breeding_date_J_snum, prop_molt_score_max, max_breeding_score, age, migratory_status) %>% distinct() %>% group_by(sex) %>% mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>% ungroup() %>% left_join(., moult_commence_ring_season %>% mutate(season = as.numeric(as.character(season))), by = c("rings_comb", "season"))```#### effect of breeding plumage score on prolonged breedingdo individuals with a fuller breeding plumage (for their sex) breed later than the rest of the population? No effect```{r}mod_max_score_last_breeding_date <- lme4::lmer(last_breeding_date_J_snum ~ std_max_breeding_score + (1| rings_comb) + (1| season),data = ind_prop_molt_scores_last_breeding_date)tbl_regression(mod_max_score_last_breeding_date, intercept =TRUE)ggplot() +geom_line(data =as.data.frame(allEffects(mod_max_score_last_breeding_date))$std_max_breeding_score,aes(x = std_max_breeding_score, y = fit)) +geom_ribbon(data =as.data.frame(allEffects(mod_max_score_last_breeding_date))$std_max_breeding_score,aes(x = std_max_breeding_score, ymax = upper, ymin = lower),lwd =1, alpha =0.25) +geom_point(data = ind_prop_molt_scores_last_breeding_date,aes(x = std_max_breeding_score, y = last_breeding_date_J_snum),alpha =1, shape =20) +theme(text =element_text(family ="Verdana"),axis.ticks =element_blank(),plot.margin =unit(c(0.2,0.2,0.2,0.2), "cm"),panel.spacing =unit(1, "cm"),panel.background =element_rect(fill ='transparent', color =NA),plot.background =element_rect(fill ='transparent', color =NA), panel.grid.major.y =element_line(size =0.25, lineend ="round", colour ="grey90"),panel.grid.minor =element_blank(),panel.grid.major.x =element_line(size =0.25, lineend ="round", colour ="grey90")) +ylab("year-scaled date of last breeding attempt") +xlab("sex-scaled max. breeding plumage score")```#### effect of breeding plumage score on early breedingdo individuals with a fuller breeding plumage (for their sex) breed earlier than the rest of the population? No effect```{r}mod_max_score_first_breeding_date <- lme4::lmer(first_breeding_date_J_snum ~ std_max_breeding_score + (1| rings_comb) + (1| season),data = ind_prop_molt_scores_first_breeding_date)tbl_regression(mod_max_score_first_breeding_date, intercept =TRUE)ggplot() +geom_line(data =as.data.frame(allEffects(mod_max_score_first_breeding_date))$std_max_breeding_score,aes(x = std_max_breeding_score, y = fit)) +geom_ribbon(data =as.data.frame(allEffects(mod_max_score_first_breeding_date))$std_max_breeding_score,aes(x = std_max_breeding_score, ymax = upper, ymin = lower),lwd =1, alpha =0.25) +geom_point(data = ind_prop_molt_scores_first_breeding_date,aes(x = std_max_breeding_score, y = first_breeding_date_J_snum),alpha =1, shape =20) +theme(text =element_text(family ="Verdana"),axis.ticks =element_blank(),plot.margin =unit(c(0.2,0.2,0.2,0.2), "cm"),panel.spacing =unit(1, "cm"),panel.background =element_rect(fill ='transparent', color =NA),plot.background =element_rect(fill ='transparent', color =NA), panel.grid.major.y =element_line(size =0.25, lineend ="round", colour ="grey90"),panel.grid.minor =element_blank(),panel.grid.major.x =element_line(size =0.25, lineend ="round", colour ="grey90")) +ylab("year-scaled date of first breeding attempt") +xlab("sex-scaled max. breeding plumage score")```#### effect of prolonged breeding on breeding plumage scoredo individuals with a late breeding attempts tend to have a fuller breeding plumage? No effect```{r}mod_last_breeding_date_max_score <- lme4::lmer(std_max_breeding_score ~ last_breeding_date_J_snum + (1| rings_comb) + (1| season),data = ind_prop_molt_scores_last_breeding_date)tbl_regression(mod_last_breeding_date_max_score, intercept =TRUE)ggplot() +geom_line(data =as.data.frame(allEffects(mod_last_breeding_date_max_score))$last_breeding_date_J_snum,aes(x = last_breeding_date_J_snum, y = fit)) +geom_ribbon(data =as.data.frame(allEffects(mod_last_breeding_date_max_score))$last_breeding_date_J_snum,aes(x = last_breeding_date_J_snum, ymax = upper, ymin = lower),lwd =1, alpha =0.25) +geom_point(data = ind_prop_molt_scores_last_breeding_date,aes(x = last_breeding_date_J_snum, y = std_max_breeding_score),alpha =1, shape =20) +theme(text =element_text(family ="Verdana"),axis.ticks =element_blank(),plot.margin =unit(c(0.2,0.2,0.2,0.2), "cm"),panel.spacing =unit(1, "cm"),panel.background =element_rect(fill ='transparent', color =NA),plot.background =element_rect(fill ='transparent', color =NA), panel.grid.major.y =element_line(size =0.25, lineend ="round", colour ="grey90"),panel.grid.minor =element_blank(),panel.grid.major.x =element_line(size =0.25, lineend ="round", colour ="grey90")) +xlab("year-scaled date of last breeding attempt") +ylab("sex-scaled max. breeding plumage score")```### Relationships between breeding schedule and moult#### effect of prolonged breeding on moultdo individuals with late breeding attempts have a delayed moult? No effect```{r}mod_last_breeding_date_prop_molt_mid <-lmer(mid_point_date_J ~ last_breeding_date_J_snum + (1| rings_comb) + (1| season),data = ind_prop_molt_scores_last_breeding_date)tbl_regression(mod_last_breeding_date_prop_molt_mid, intercept =TRUE)ggplot() +geom_line(data =as.data.frame(allEffects(mod_last_breeding_date_prop_molt_mid))$last_breeding_date_J_snum,aes(x = last_breeding_date_J_snum, y = fit)) +geom_ribbon(data =as.data.frame(allEffects(mod_last_breeding_date_prop_molt_mid))$last_breeding_date_J_snum,aes(x = last_breeding_date_J_snum, ymax = upper, ymin = lower),lwd =1, alpha =0.25) +geom_point(data = ind_prop_molt_scores_last_breeding_date,aes(x = last_breeding_date_J_snum, y = mid_point_date_J),alpha =1, shape =20) +theme(text =element_text(family ="Verdana"),axis.ticks =element_blank(),plot.margin =unit(c(0.2,0.2,0.2,0.2), "cm"),panel.spacing =unit(1, "cm"),panel.background =element_rect(fill ='transparent', color =NA),plot.background =element_rect(fill ='transparent', color =NA), panel.grid.major.y =element_line(size =0.25, lineend ="round", colour ="grey90"),panel.grid.minor =element_blank(),panel.grid.major.x =element_line(size =0.25, lineend ="round", colour ="grey90")) +xlab("year-scaled date of last breeding attempt") +ylab("mid-point date between\nlast max score photo and first photo during moult")```#### effect of prolonged breeding on prolonged retention of full plumagedo individuals with late breeding attempts retain their full breeding plumage longer?There is a positive relationship between the date when an individual was last seen with its maximum breeding plumage and the date when it was last observed breeding.```{r}mod_last_breeding_date_prop_molt_latest <-lmer(latest_max_date_J ~ last_breeding_date_J_snum + (1| season),data = ind_prop_molt_scores_last_breeding_date)tbl_regression(mod_last_breeding_date_prop_molt_latest, intercept =TRUE)ggplot() +geom_line(data =as.data.frame(allEffects(mod_last_breeding_date_prop_molt_latest))$last_breeding_date_J_snum,aes(x = last_breeding_date_J_snum, y = fit)) +geom_ribbon(data =as.data.frame(allEffects(mod_last_breeding_date_prop_molt_latest))$last_breeding_date_J_snum,aes(x = last_breeding_date_J_snum, ymax = upper, ymin = lower),lwd =1, alpha =0.25) +geom_point(data = ind_prop_molt_scores_last_breeding_date,aes(x = last_breeding_date_J_snum, y = latest_max_date_J),alpha =1, shape =20) +theme(text =element_text(family ="Verdana"),axis.ticks =element_blank(),plot.margin =unit(c(0.2,0.2,0.2,0.2), "cm"),panel.spacing =unit(1, "cm"),panel.background =element_rect(fill ='transparent', color =NA),plot.background =element_rect(fill ='transparent', color =NA), panel.grid.major.y =element_line(size =0.25, lineend ="round", colour ="grey90"),panel.grid.minor =element_blank(),panel.grid.major.x =element_line(size =0.25, lineend ="round", colour ="grey90")) +xlab("year-scaled date of last breeding attempt") +ylab("last date photographed with max score")```### Relationships between age and breeding phenology#### effect of age on prolonged breeding do older individuals conclude their breeding season later? No```{r}mod_age_last_breeding_date <- lme4::lmer(last_breeding_date_J_snum ~ age + (1| rings_comb) + (1| season),data = ind_prop_molt_scores_last_breeding_date)tbl_regression(mod_age_last_breeding_date, intercept =TRUE, label =list(age ~"Age"))ggplot() +geom_line(data =as.data.frame(allEffects(mod_age_last_breeding_date))$age,aes(x = age, y = fit)) +geom_ribbon(data =as.data.frame(allEffects(mod_age_last_breeding_date))$age,aes(x = age, ymax = upper, ymin = lower),lwd =1, alpha =0.25) +geom_point(data = ind_prop_molt_scores_last_breeding_date,aes(x = age, y = last_breeding_date_J_snum),alpha =1, shape =20) +theme(text =element_text(family ="Verdana"),axis.ticks =element_blank(),plot.margin =unit(c(0.2,0.2,0.2,0.2), "cm"),panel.spacing =unit(1, "cm"),panel.background =element_rect(fill ='transparent', color =NA),plot.background =element_rect(fill ='transparent', color =NA), panel.grid.major.y =element_line(size =0.25, lineend ="round", colour ="grey90"),panel.grid.minor =element_blank(),panel.grid.major.x =element_line(size =0.25, lineend ="round", colour ="grey90")) +xlab("age") +ylab("year-scaled date of last breeding attempt")```#### effect of age on early breeding do older individuals start breeding earlier? No```{r}mod_age_first_breeding_date <- lme4::lmer(first_breeding_date_J_snum ~ age + (1| rings_comb) + (1| season),data = ind_prop_molt_scores_first_breeding_date)tbl_regression(mod_age_first_breeding_date, intercept =TRUE, label =list(age ~"Age"))ggplot() +geom_line(data =as.data.frame(allEffects(mod_age_first_breeding_date))$age,aes(x = age, y = fit)) +geom_ribbon(data =as.data.frame(allEffects(mod_age_first_breeding_date))$age,aes(x = age, ymax = upper, ymin = lower),lwd =1, alpha =0.25) +geom_point(data = ind_prop_molt_scores_first_breeding_date,aes(x = age, y = first_breeding_date_J_snum),alpha =1, shape =20) +theme(text =element_text(family ="Verdana"),axis.ticks =element_blank(),plot.margin =unit(c(0.2,0.2,0.2,0.2), "cm"),panel.spacing =unit(1, "cm"),panel.background =element_rect(fill ='transparent', color =NA),plot.background =element_rect(fill ='transparent', color =NA), panel.grid.major.y =element_line(size =0.25, lineend ="round", colour ="grey90"),panel.grid.minor =element_blank(),panel.grid.major.x =element_line(size =0.25, lineend ="round", colour ="grey90")) +xlab("age") +ylab("year-scaled date of first breeding attempt")```## Meeting notes:### 13-Nov-24: Bart, Luke, and Bashar via zoom- add the breeding data to the analysis: for each individual include the date of last breeding (i.e., last date seen with nest or brood) as a predictor of molt timing. Also do an analysis relating the timing of last breeding to the age of the individual. Bashar will send Luke the breeding data he received from Ted.- drop pre-October data from the molt timing analysis.- next steps: Bashar to eventually write up this study as a report (which could eventually be used as his MSc thesis if he wants). Priority will be on him writing up a proposal for his thesis (which has a concrete deadline and is graded). Bashar will present the study at our weekly Monday meeting in the first few months of 2025.